INTERACTIVE COMPUTATION OF RADIATION VIEW FACTORS* 


A.F. EMERY 
H.R. MORTAZAVI 
C.J. KIPPENHAN 

UNIVERSITY OF WASHINGTON 
SEATTLE, WASHINGTON 


ABSTRACT 

The development of a pair of computer programs to calculate the radiation 
exchange view factors is described. The surface generation program is based 
upon current graphics capabilites and includes special provisions which are 
unique to the radiation problem. The calcul ational program uses a combination 
of contour and double area integration to permit consideration of radiation 
with obstructing surfaces. Examples of the surface generation and the 
calculation are given. 


INTRODUCTION 

The calculation of the radiation exchange between two surfaces by the 
usual engineering method 
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first requires the evaluation of the view factor, F-|2 » which represents the 
fraction of energy leaving one diffuse surface which is intercepted by the 
second surface. F^2 is defined by[l] 
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Although the equation is simple in appearance, its evaluation is fraught 
with difficulties since it requires a full and precise description of both 
surfaces. In addition, if the view from surface 1 to surface 2 is obscured by 
an interposed surface or object, some provisions must be made to account for 
that porion of the view which is occluded. It js important to recognize that 
the obstructed view is not a constant, but varies, depending upon the position 
of the elemental area dA-| . 

The calculation of the view factor is thus really two problems in one: 


1. To devise an efficient way of describing the surfaces and their 
spatial relationship to each other. 

2. To calculate the resulting view factors. 


This paper describes the development of a pair of computer programs to do 
this. The calculational program was developed first and used with hand input 
in a batch mode. It later became apparent that efficient use of the program 
required some form of automatic surface generation and the appropriate 
subroutines were added. It also became painfully obvious that the average user 
made so many errors (both in key punching and in defining the surfaces) that an 
interactive program, with graphic capability, was necessary. 

SURFACE GENERATION 

The first impression was that a standard CAD/CAM program could be easily 
modified to provide the needed interactive capability. However, some of the 
unique requirements of the view factor calculation program led us to conclude 
that a special program would be more appropriate. This program was based upon 
the following assumptions: 

1. Surfaces would be only 3 or 4 sided with straight sides 

2. To simplify the obstruction calculations, all surfaces would be planar 

3. A surface would radiate from one side only. Thus every plate, no 
matter how thin, would require 2 surfaces. 

4. The direction of the surface normal must be uniquely and simply 
defined. 

5. Surfaces could not penetrate one another. 


Some of these requirements were needed to calculate the view factors. 
Others were imposed to simplify the surface generation. Several additional 
requirements were found to be useful, although not absolutely necessary, 
namely: 
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1. All surfaces are initially generated in the x,y plane, facing upwards. 
Three dimensional surfaces (cylinders etc) would be oriented along the 
z axis 

2. Only a global coordinate system would be used. Individual surfaces 
would not carry along an embedded coordinate system, 

3. Surface corner nodes would always be numbered in a counter clockwise 
direction to specify the surface normal. 

4. Curved surfaces would be represented by a combination of triangles and 
quadrilaterals, the size of each adjusted by the user to give the best 
representation. 


The different capabilities of the generation program are illustrated in 
figure 1. Because radiation is a surface phenomena, the surface orientation 
(i.e., the direction of the surface normal) must be uniquely defined. In 
testing the program, particularly with inexperienced users, we have found that 
it must not only be interactive but must permit the user to modify any command 
at any tTime and must have a very complete graphics capability. The following 
sections describe the different features of the program. 

A. Surface Types 

Three different surface types are used: minor, major and groups. A major 
surface is one that is generated by a single command. Typical major surfaces 
are listed in table 1. A major surface can be manipulated as a single entity. 
A group is a collection of major surfaces and is also a single entity. When 
the major surface is generated, it may be composed of many minor or 
sub-surfaces. Figure 2 illustrates the development of a group and shows the 
minor surfaces. Minor surfaces cannot be treated as independent entities since 
this would lead to possible erratic and unacceptable distortions of the major 
surfaces. 


TABLE 1 

Typical Major Surfaces 


Two Dimensional 

1. Triangle 

2. Rectangle 

3. Quadrilateral 

4. Circle (annulus, sector) 

5. Ellipse (annulus, sector, orientation) 


Three Dimensional (inside or outside radiation, top and/or bottom) 

1. Box 

2. Cylinder (right, slant, annular, sector) 

3. Cone (right, slant, frustrum, annular, sector) 

4. Sphere (sector, cutting planes) 
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B. Command Data Base 


Because even experienced users tend to make many mistakes or desire to 
make a substantial number of modifications, it is important that all commands 
be stored for future use and modification. Commands are either informative or 
functional. All functional commands are stored in a data base and the user 
can: 

1. List any or all of the commands 

2. Insert new commands between existing commands. 

3. Delete commands 

4. Temporarily suspend commands 

5. Repeat the commands 

6. Interrupt the repeat 

7. Call for information or graphical display at any time during the 
repetition of the commands. 

Since surface movements are non-commutative and because the current structure 
cannot be generated by only a portion of the commands, any changes in the data 
base must be accomplished by restarting the command sequence (i.e. repeating) 
and modifying during the subsequent generation process. 


C. Surface Manipulation:; 

The user must be able to manipulate each entity (major or group) by such 
movements as: 

1. Translation 

2. Rotation 

3. Replication (i.e, duplicating an original set of surfaces and 
manipulating the new set) 

4. Scaling in x,y or z coordinates 


In addition, all surfaces (including minor surfaces) must be subject to: 

1. Deleting (temporary or permanent) 

2. Adding 

3. Restoring (if previously deleted) 

4. Joining to another surface to form a new entity 

5. Separating from another 

6. Reversing the normal direction 

7. Combining with others to form a single radiating surface but still 
capable of independent manipulation. 
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D. Graphical Display 

The key to effective surface generation is a high speed graphical display. 
Our experience has been that the typical user will call for a display after 
every one or two functional commands. Furthermore, the user will generally ask 
for more than one view. Thus the graphics must be fast and versatile. The 
minimum number of graphical views and commands appears to be: 

1. Orthographic views (top, bottom, left, right, front, back) singly or in 
groups 

2. Perspective or isometric views 

3. Hidden line removal from any or all views 

4. Variable field of view 

5. Rotation and translation of the object. 

6. Deletion of selected surfaces from the view 

7. Numbering of selected surfaces or corner nodes 

8. Display of surface normals 

9. Views as seen from one surface to another 


Although the isometric view gives the best overall picture, the 
orthographic views prove to be the most useful when precise movements are 
necessary. Because obstructions so completely change the radiation exchange, 
it is critical that the user be able to look from any one surface to any other 
to check for obstructions. This must be done with hidden line removal to yield 
the maximum information. 

Figure 3 illustrates the incorrect movement of panels on a cylinder as 
shown by an isometric view and an orthographic view — the latter being necessary 
for the correct placement. 

Although a wire frame display is useful, the usual structure has so many 
surfaces that a hidden line removal is necessary. Unfortunately hidden views 
are very time consuming and consequently the program incorporates a somewhat 
inexact, but very fast, routine. Part of its speed is based upon the use of 
planar, straight sided surfaces with counter clockwise numbering. Any other 
surface requires some type of raster scan hidden surface algorithm with an 
unacceptable increase in computing time[2]. It is interesting to note that a 
fast graphics display may not give the desired information. We have found that 
if the lines are drawn slowly enough so that the viewer can get a feeling for 
the order of generation, it is easier to visualize the structure. On the other 
hand, if the display is nearly instantaneous, many viewers cannot recognize the 
structure. This is particularly true if an old command data base is being 
reused. Finally the surface normals must be displayed since a surface has only 
one active side. Most user errors appear to be related to an incorrect 
orientation of a surface after a series of manipulations have been made. 
Figure 4 illustrates a set of double surfaces to represent .panel s and the 
erroneous orientation of the replicated set is clearly seen. 
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E. Surface Refinement 


After the basic structure has been constructed, three additional functions 
are needed: 

1. Refinement by subdividing surfaces 

2 . Assuring planarity of surfaces 

3. Condensing nodes 


Refinement is simply the subdivision of surfaces into triangular 
sub-surfaces. Since triangles are always plane, planarity is assured by 
dividing any non-planar quadrilateral into 2 triangles. Obviously, refinement 
also guarantees planarity. Finally, upon exit any two nodes which are within 
a prescribed distance e are condensed into one node*, all inactive surfaces and 
unused nodes are deleted. The data file is created for use by the view factor 
calculating progam and for subsequent plotting. The command data base is 
closed for future use. 


CALCULATING THE VIEW FACTORS 

Because analytical expressions for F- ^ exist only for simple surfaces [3], 
the computation of the view factors for complex structures is practical only by: 

1. Double area integration 

2. Contour integration 

3. Nusselt projection 

4. Monte Carlo 


The Monte Carlo method[4] tends to be too expensive for most surfaces even 
if one of the adaptive techniques [5] is used. Nusselt projection (i.e., 
the use of the unit sphere) generally produces a curved projected surface and 
consequently requires a complex integration to find the projected area. In 
addition, obstructions will distort the contour of the viewed surface, 
rendering the calculation of the projected edges very difficult. In this case, 
one simply projects the position of an elemental area, not the edges of the 
contour. Since the projection must be done from every point on the viewing 
surface, it essentially reduces to the double area integration method. 
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A. Double Area Integration 


The double area integration is performed in two steps (figure 5), 
1. Mapping the surface onto a unit square by 


~ I Ni (?,n) 


( 3 ) 


The functions N^- are the usual finite element isoparametric 
functions [6]. 

2. Evaluating the integral 



where (?,n) are the Gaussian points in the unit square, are 
the Jacobian of transformation evaluated at these points, nj^ are the 
Gaussian weights and the superscripts i and j refer to the surfaces 
and 
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B. Contour Integration 

The double area method is the best of the three, but it is not 
sufficiently accurate when portions of the surfaces are close to each other 
since the denominator of equation 1 becomes very small and the integral nearly 
singular. For these reasons we adopted the contour integration method as 
described by Mitalas and Stephenson [7]. In this method, the double area 
integration can be replaced by 
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where the various terms are defined on figure 6 and 

D(k,£) = + n.n. {I m n = direction cosines) 

This formulation is not sensitive to the separation distance of the edges since 
the term is not singular as the distance r goes to zero. We note that 

planar or non-planar surfaces can be treated equally well. 

C. Comparison 

In general , when numerically implemented, either Contour integration or 
Double Area integration will produce acceptable results. However, when two 
surfaces have a common edge (the adjoint problem) the Double Area integration 
method may perform very poorly[8, 9]. Table 2 lists comparable values obtained 
by the two methods for the surface shown on figure 7. Because the greatest 
portion of the radiation occurs in the corner, where the surfaces are the most 
proximate, the Double Area method tends to be inaccurate. By contrast the 
Contour integration method is very accurate, even with very few edge elements. 
When the surfaces are separated slightly, there is some improvement in the 
Double Area results, but not sufficient to justify its choice over the contour 
method. 


TABLE 2 

Percentage Error in the Numerical Calculation of the View Factor 
between Two Surfaces of Equal Breadth (L=H) 

(see figure 7) 


Contour Integration Double Area Integral 

Infinite Finite Area Infinite Finite Area 

Strip Strip 


d/L 

S=0.0 

S=0.0 

S=0.0 

S=0.1 L 

S=0.0 

S=0.1 L 

1,0 

0.2% 

3.7 

20.7 

21.5 

59.1 

57.6 

0.5 

0. 

0.8 

12.9 

11.2 

32.0 

24.1 

0.3 


0.4 

9.1 

6.6 

21.8 

12.1 

0.2 


0.1 

5.7 

2.8 

13.2 

4.9 

0.1 


0. 

2.9 

0.7 

6.7 

1.2 

0.05 



1.5 

0.2 

3.4 

0.3 


d = 

dx=dy=dz for the 

Double Area 

Integration 




edge segment length for the Contour Integration 
S = Separation distance along the x coordinate 


We note that the use of the Contour integral with as few as 5 segments on an 
edge gives acceptable results, while the use of the Double Integral with a 
corresponding spacing is quite unacceptable. 
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D. Program Structure 


The basic program structure is illustrated in Table 3, 


TABLE 3 
PROGRAM FLOW 

1. Geometrical Input (definition of surfaces) 

2. Unobstructed View Factor Computations using contour integration 

3. Obstructed View Factor Computations 

1. determining if two surfaces, i and j, can see each other 

2. determining if a 3rd surface is interposed between the pair 

3. subdividing surfaces i and j into elemental areas 

4. constructing a ray between points on surfaces i and j and 
determining if it is intercepted 

5. computing F.. using the double area integral 


Each of the three major functions of the program will be described in 
detail in the following sections. 


D1 . Input 

Data input is from the output file of the Surface Generation program. 


D2. Unobstructed View Factor Calculations 

The values of Fy are computed directly by using the Contour Integration 
equation. Calculational time is reduced by the use of the reciprocity 
relationship, at the cost of the loss of some additional information. 

D3. Obstructed View Factor Calculations 

(A) Elimination of non-obstructing Surfaces 

The calculation of F-rj when one or more surfaces are interposed between 
surfaces i and j is tne most difficult and time consuming part of the 
calculation. In order to accelerate this calculation, every possible effort 
must be made to eliminate all non-obstructing surfaces from consideration. 
Consider the situation schematically shown on figure 8. The best calculational 
procedure found to date is: 
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I I II II 


II III II 


1. Using the original x,y,z coordinate system, define the smallest 

rectangular parallelepiped (the view prism) which contains surfaces i 
and j (indicated by the dashed lines). Eliminate all surfaces which 
do not penetrate this view prism. This process is usually referred to 
as "clipping". Check to see if any surface completely fills the view 
prism, since if the view is totally obscured no further computations 
are needed. _ _ 

2. Define a new coordinate system, 'x,y,Y, with the T axis directed along 

the line connecting the centroids of surfaces i and j and eliminate 
all possible surfaces by another clipping pass. Again test for total 
blockage. _ _ 

3. Define a third coordinate system, x,y by rotating the x,y coordinate 
system until the rectangular area which encloses both surfaces i and j 
is a minimum (figure 8b) and perform another clipping test and test 
for total blockage. 


It cannot be emphasized too strongly that acceptance or rejection of a 
possible obstructing surface can only be done by the siinple clipping operation 
associated with view prisms which are constructed from rectangular 
parallelepipeds. Any other procedure to test for penetration of the view prism 
calls for geometric calculations which are unacceptably expensive. Since the 
transformations involved in steps 2 and 3 require matrix operations on all of 
the corner coordinates of all candidate obstructions, it is imperative that 
each successive step eliminate as many obstructing surfaces as possible. We 
have considered using the graphic capability of the generating program to 
detect the obstructing surfaces, but it has proven to be too time consuming and 
i nefficient. 

Once the final set of possible obstructing surfaces has been determined, 
F^; can be calculated in two ways. Consider the view of surface j as seen from 
a segment on the contour of surface i as shown in figure 9 . The obstructing 
surfaces obscure part of surface j, leading to the formation of the two visible 
sub-areas indicated by the dotted line contours. We note that: a) the 
determination of the unobstructed part of the surface is the classical hidden 
surface problem of graphic display for which, currently, there are no efficient 
methods; b) the number of non-contiguous visible sub-areas and the number of 
straight line segments encompassing these sub-areas may vary considerably as 
seen from different points of the contour of surface i. For this reason. 
Contour Integration is not an acceptable method; hence, we must utilize Double 
Area Integration. 
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(B) Ray Interception 
We express as 




"i "j 

V I(l<,£)dA,^dA^ 


( 8 ) 


where f|^£= 1 if the ray between dA|^and dAg^is unobstructed 
= 0 if obstructed. 


Since fj^^is a discontinous function of position, the higher accuracy of 
Gaussian .quadrature is not always realized; and since Gaussian quadrature 
requires substantially more numerical operations, it has proven best to use 
Newton-Cotes integration. 

Once each surface has been divided into elemental areas, the centroids of 
these areas determined, and the rays between each of the points on surface i to 
each of the points on surface j defined, the calculation of proceeds as 

follows (figure 10); 

1. Determine if the angle between the outward normals to surfaces i and j 

and the ray Piq^is greater than 90 degrees. If so, set fk£=0 since the 

surfaces cannot see each other. Note that this must be done for every 
ray, since the inability of one point on surface i to see any given 
point on surface j does not ensure that other parts of the two 

surfaces cannot see each other. 

2. Determine the intersection of the ray Pj^^with each of the possible 
obstructing surfaces by examining each in turn. This is done by: 

1. Finding the intersection of the ray Pkg,with the infinite surface 
which contains the surface. Because the intersection requires an 
iteration for arbitrary surface, but not for planes, only plane 
surfaces are permitted. The determination of the intersection is 
best effected by pre-calculating the equation of the plane of each 
surface and the transformation R 


fx"' 


fx] 

V 

U^J 

■ = [R] ’ 

A 


(9) 


which produces the coordinate system x',y',z' for which x' and y' 
are in the plane of the obstructing surface and z' normal to it. 
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2. Determining if the point of intersection is within the surface or 
not. If it is, if Determining whether the 

intersection point is within the surface by mapping the surface to 
a square and checking the values of 5 and n is not efficient 
because the non-linear mapping requires an iterative solution for 
5 and n . The most efficient way appears to be by computing the 
angle between the line drawn from the corner node to the 
intersection point and between the corner node and the next corner 
node. For convex surfaces, if this angle is negative for any 
corner, the point is not within the surface and the ray is not 
blocked. 


The most efficient calculation is one in which the obstruction is detected 
early, thus eliminating the consideration of the remaining surfaces. This is 
only possible if the candidate obstructing surfaces can be sorted in the order 
of their included angle as seen from the elemental area, dA|,. Unfortunately, 
this determination of the included angle is equivalent to finding the view 
factor from the elemental area dA|^ to that portion of the obstructing surface 
which is within the view prism between points dA|. and dAo . Thus far we have 
not found a sorting method which gives a net reduction in computational times. 
Our experience is that maintaining a fixed order of checking the rays for 
interception in a single view prism is an effective method. (It should be 
noted that this entire problem can be compared to the usual graphical display 
problem for perspective views of hidden surfaces, but since in essence the 
viewing point must range over the entire surface i, calculational ly it is the 
square of this classical problem.) 

That the Double Area Integration is inaccurate for surfaces with a common 
edge has already been discussed; yet only the Double Area Integration is useful 
for the obstructed view calculations. Therefore it is important to avoid this 
situation by appropriately defining the surfaces. 


EXAMPLES 

We present two examples of the use of the program. Figure H represents 
an enclosure with a dividing panel, with the surfaces as numbered. The 
dividing partition must be expressed as two surfaces, infini tesimally separated, 
not one. Because surface 1 has an obstructed view of surface 9 with which it 
has a common edge, it must be represented as two surfaces, lA and IB. In this 
way, surface lA which adjoins surface 9 has an unobstructed view and is treated 
by Contour Integration. Surface IB, which has an obstructed view, is not 
adjoining and can be successfully treated by Double Area Integration. Figure 
12 represents an enclosed cylinder which obstructs the view of the other 
surfaces. Table 4 lists typical execution times. The obstruction calculations 
were carried out using the three different elimination methods and different 
numbers of rays per surface and the accelerating effect of the clipping 
routines is clearly shown. Because of the simple geometries and the 
orientation of the surfaces along coordinate axes, the second clipping was not 
effective since the slight reduction in computational time due to the 
elimination of a few additional surfaces was less than the time necessary to 
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perform the coordinate transformations. When the x,y axes are rotated by about 
22 degrees, the first clipping effectiveness is reduced by about 30% and the 
value of the second clipping is more pronounced. For general problems, all 
clipping procedures must be used. 


TABLE 4 

Calculation Times 

for the Example Problems depicted in Figures 11 and 12 


Adjoint Surfaces 
figure 10 
(12 surfaces) 


Enclosed Cylinder 
figure 11 
(48 surfaces) 


Number of rays 

9 16 

9 

16 

16 

per surface 




(rotated 

coordinates) 

Cl ipping 

Calculational 

Times 

in Central Processor Seconds 

routine 

on CDC Cyber 

175/750 


none 

5.6 14.8 

99.0 

274.1 

274.1 

1st pass 

2.9 6.3 

35.7 

82.8 

110.8 

2nd pass 

3.0 6.3 

38.9 

82.5 

82.5 

2nd pass(opt=2) 




76.3 


CONCLUSIONS 

Our experience in using the programs has emphasized that radiation 
view factors of complex structures can only be accomplished if a fast, 
interactive program with good graphics capability is used. Although it would 
be more satisfying if curved surfaces could be treated, high program efficiency 
requires that they be modelled by an assemblage of flat triangles or 
quadril aterals. 


The view-factor calculational algorithm described has proven to be very 
effective for surfaces which have an obstructed view of each other. Under some 
conditions, statistical methods may prove to be more efficient, but for general 
configurations the combination of the Contour Integration and Double Area 
Integration methods, in concert with ray intersection calculations, has proven 
to be about the fastest method currently available. Further increases in 
computing efficiencies are possible if hardware perspective view, hidden 
surface devices are used, but such devices are not currently available for 
digital computers used in thermal analyses. From another point’ of view, the 
use of the rays may be regarded as a highly adaptive form of the Monte Carlo 
method which bears the same relationship to the usual method that the 
quasi-determini Stic Exodus [10] method bears to the usual Monte Carlo method 
for solving diffusion problems. 
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NOMENCLATURE 

A-j Area of surface i 

dAj Elemental area of surface 

F-jj Vi ewf actor from surface i to surface j 

hjj Gaussian weight for point k in surface i 

J(P|^) Jacobian of transformation at point k 
n,- Number of elemental areas in surface j 

fT Isoparametric shape functions 

P|^ Gaussian point k on surface i 

r-jj Distance from surface i to surface j 
R Rotation matrix 

x,y,z Global coordinate system 

x.y.T Rotated coordinate system 

0 Angle between r^j and a surface normal 

Surface coordinates of the unit square. 
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Figure 1.- Surface generation-program capabilities. 




Figure 2.- Development of a group of surfaces 






Figure 3. 


Fi gure 




Incorrect placement of panels, isometric and orthographic views. 



Display of surface normals showing incorrect orientation. 
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Figure 5.- Double area integration (showing mapping onto a unit square). 




Figure 6.- Definition of angles and lengths for contour integration. 



I 



Figure 7.- View factors between surfaces at right angles to each other and separated 

by the distance S. 
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Figure 9.- View of surface i with interposed obstructions showing the multiple 

segmented contours. 


INFINITE PLANE 
CONTAINING OBSTRUCTION 



Figure 10.- Geometry for determining if the ray P(k,l) intercepts an obstructing 

surface. 
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